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ABSTRACT 

Eighty planetary systems of two or more planets are known to orbit stars other than the Sun. For 
most, the data can be sufhciently explained by non-interacting Keplerian orbits, so the dynamical 
interactions of these systems have not been observed. Here we present 4 sets of lightcurves from 
the Kepler spacecraft, which each show multiple planets transiting the same star. Departure of the 
timing of these transits from strict periodicity indicates the planets are perturbing each other: the 
observed timing variations match the forcing frequency of the other planet. This confirms that these 
objects are in the same system. Next we limit their masses to the planetary regime by requiring the 
system remain stable for astronomical timescales. Finally, we report dynamical fits to the transit 
times, yielding possible values for the planets' masses and eccentricities. As the timespan of timing 
data increases, dynamical fits may allow detailed constraints on the systems' architectures, even in 
cases for which high-precision Doppler follow-up is impractical. 

Subject headings: planetary systems; stars: individual (KID 10358759 / KOI-738 / Kepler-29; KID 
3832474 / KOL806 / Kepler-30; KID 9347899 / KOL935 / Kepler-31, KID 9787239 
/ KOI-952 / Kepler-32); planets and satellites: detection, dynamical evolution and 
stability; methods: statistical 
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1. INTRODUCTION 

So far, 170 systems with more than one transit- 
ing planet candidat e have been discovered by Kepler 
(|Borucki et al.|[2i311[ ). enabhng simple dynamical mod- 
els for an astoundin g variety of planetary systems 
(|Lissauer et alj l2011bD . In particular, two very im- 
portant dynamical quantities, the planetary period and 
phase, are measured to high precision by Kepler. If we 
further assume, as is true for the Solar System, that the 
planets orbit in the same direction, in nearly the same 
plane, and on nearly circular orbits, then we have a fidu- 
cial model that specifies the positions of the planets as 
a function of time. However, planets do not follow inde- 
pendent Keplerian orbits; instead planets interact with 
each other gravitationally. These small perturbations to 
the orbits affect the transit times, which can be calcu- 
lated via numerical simulations. If the calculated transit 
timing variations match essential aspects of the observed 
data, then we obtain independent confirmation that two 
objects are orbiting the same star. 

In this paper, we develop simple physical models 
and compare them to Kepler data to confirm 4 plan- 
etary systems, Kepler-29-32. Instead of exhaustively 
modeling t he systems, a. s has been our prac t ice up 
until now (iHolman et al.l 120101: iLissauer et all l2011al: 
ICochran et al.ll2niin . we content ourselves with provid- 
ing upper limits to the planetary masses, based on the 
principle that the system is extremely unlikely to be dy- 
namically unstable on timescales much less than the age 
of the star. Thus, we infer planetary masses for 9 ob- 
jects in these 4 systems, and deem them planets on this 
basis. This paper is being published contemporaneously 
with two companion papers, which explore anticorrelated 
transit timi ng variations for t he confirmation of pairs of 
planets, bv iFord et al.l ()2012D (hereafter Paper II) and 
iStefFen et alT i 20121 ) (hereafter Paper III) , which together 
with this paper form a catalog of 10 confirmed systems 
of multiple-transiting planets. 

The organization of this paper is as follows. Sec- 
tion [2] introduces the target stars, describes the transit 
timing data, and discusses how we assess whether the 
transit timing signal varies on the theoretically-expected 
timescale. Section 13.11 applies that methodology to con- 
firm 4 planetary systems, and section 13.21 discusses the 
identification of their host stars. Section |4] is devoted to 
constraining the masses of the planets, primarily by dy- 
namical stability, but also via modeling the transit timing 
signal itself. Section [5] closes with a discussion of these re- 
sults and a perspective about the future of transit-timing 
confirmation of multiply-transiting systems with Kepler. 

2. METHODS 
2.1. Determination of Stellar Parameters 

Identifying properties of the host stars studied here are 
in Table [U e.g., their names in various catalogs. 

For stars Kepler-29, Kepler-30, and Kepler-32, we ob- 
tained spectra from several different telescopes as part 
of the Kepler team's regular ground-based follow-up 
program. Two independent analyses were performed, 
a standard analysis using 'Spectroscopy Made E asy' 
(jValenti fc Piskunovl [1991 iValenti fc Fischeii [20051) . as 
well as a newly-formulated analysis called SPC (Buch- 
have et al. in prep.), to obtain the parameters Teg, logg. 



[Fe/H], and when possible wsini. For Kepler-31 we did 
not obtain a spectrum, but instead a dopted these values 
from the Kepler input catalog (KIC; [Brown "eFaI1[20Tll) 
- i.e., based on color photometry - with generous error 
bars. From these values we performed a Bayesian anal- 
ysis using Yonsei-Yale stellar isochrones LYi et al.ii200l[ ) 
to extract b est-fit and 68% confidence regions on 
and (e.g. iPont fc Eve3[200l [Takeda et al.|[2007[) . All 
these values are reported in Table [2] 

2.2. Measurement of Transit Times and Parameters 

The data we use for this study are long-cadence 
lightcurves from quarters 1 through 6 (and additionally 
for Kepler-31 only, quarters 7 and 8), available at the 
Muhimission Archive at STScI (MASlE3). 

For most of the sample, we used timing results from the 
ge neral-purpose rou tines of J. Rowe, previously described 
in lFord et afl (j201lD [Paper I]. For Kepler-30b, extremely 
large variation in the transit times caused the general- 
purpose algorithm to fail (a linear model for the times of 
transit incorrectly predicts some observed times by over 
half a day, see section [3".1.2|) . Thus, arou nd each transit of 
this p lanet, we fit model transit curves (jMandel fc Agoll 
I2OO2I) and minimized the sum-of-squared residuals, the 
standard function. We sampled over a grid of 
timing offsets, and fit a parabola to the region within 
Ax^ < 7 of the minimum value. That parabola thus has 
a width which is more stable to noise properties than 
the local curvature is, and we adopt its minimum as the 
transit time and its width (where Ax^ = 1) as the error 
bar. 

The transit times for all the planets analyzed herein 
are given in Table [S] 

We also fit transit curves for all our candidates to de- 
termine best-fit parameters and errors. The first step 
was to detrend the lightcurve; we masked out each tran- 
sit and fit a polynomial over timescales of 1000 min- 
utes and divided the whole lightcurve by it, to obtain 
a detrended, normalized lightcurve /. Second, we de- 
fined a contamination-corrected lightcurve via /coir = 
(/ — c)/(l — c), where c is the fractional amount of light 
leaking into the aperture from known nearby stars, the 
"contamination" reported on MAST for each target for 
each quarter. Third, we assembled a phased lightcurve 
by subtracting the previously-measured transit time from 
each transiQ. We used only non-overlapping transits, 
for which another planet's transit center was not within 
5 hourf[3- To construct a fiux value of the long cadence 
(29.4 min), we took 20 samples e venly spaced over that 
timespan (using OCCULTSMALL; IMandel fc A"^ [200I) 
and averaged them. Finally, we fit the following four pa- 
rameters: radius ratio Rp/Ri,, duration Tdur, scaled im- 
pact parameter b, and linear limb-darkening coefficient u. 
The duration is defined as the time between mid-ingress 
and mid-egress, when the center of the planet passes over 

'^^ http:/ /archive. stsci.edu/kepler/ 

22 For the very small candidates KOI-935.04 and KOI-952.05, 
individual transit times were not reliably measured (hence they 
are not included in Table |6ll , so we used a linear ephemeris for 
fitting their transit parameters, and held their impact parameters 
b fixed at and their limb-darkening coefficients u fixed at 0.5. 

■^^ We also discarded Kepler-31c's transits near t — 245000 = 
16.8 and 400.5, as these had corrupted data, apparently due to 
detrending problems near gaps. 
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the limb of the star. We report these parameters in Ta- 
ble H 

As a check on these stellar and planetary parameters, 
and as a first indication that all the planet candidates 
in each system are transiting the target stars, we com- 
puted for each planet candidate the transit duration that 
would be obtained by an orbit that is edge-on {b = 0) 
and circular (e = 0). These values are plotted against 
the measured values in Figure [TJ We see only a little 
evidence that these computed values are slightly larger 
than the measured ones, an indication that either (i) the 
computed stellar radii are too large or stellar masses are 
too small, or (ii) the planets' orbits are systematically 
seen near pericenter of moderately eccentric orbits, or — 
most probably — (iii) the planets' orbits are not quite 
edge-on (6^0). 

2.3. Dynamical Expectations and their Observational 
Confirmation 

We begin N-body integrations near the center of the 
dataset ( t [BJD] = 2455220 ), using coplanar, circular or- 
bits that match the observed periods and transit epochs, 
and using nominal masses Mp = Mm(i?n /i?m)^ °^. This 
set of p hysical models was introduced bv lLissauer et al.l 
(|2011bl) . and it helps us diagnose the most prominent 
physical interactions. For systems of more than two plan- 
ets, we break the problem into planet pairs, in an effort 
to probe the individual inter actions. Fi r st, we calculate 
transit times as described by iFabrvckvl (|201Cl( ) from the 
beginning of the dataset until 10 years later. Next, the 
deviations of those simulated transit times from a lin- 
ear ephemeris calculated on those times (i.e., the Simu- 
lated minus Calculated diagram, S ~ C) is Fourier an- 
alyzed to find the dominant frequency (jScargld 119821 : 
iZechmeister fc Kur ster 2009), which are given in TablelH 

We test for variations at this frequency in the observed 
timing (i.e., the Observed minus Calculated diagram, 
O — C). The frequency is already determined theoret- 
ically (above), but we allow amplitude and phase of the 
sinusoid and a constant offset to vary freely to fit the 
data. This approach contrasts to a periodogram search, 
which tests many frequencies. By fixing this important 
parameter rather than scanning over it, our method is 
more sensitive than a periodogram search at finding vari- 
ations at the predicted frequency. The reason we treat 
the phase as a free parameter is that it is expected to have 
a complex dependence on the eccentricities of both the 
perturbing planet and the perturbed planet. That is, we 
cannot rely on the values from a model which assumes 
circular orbits. The same is true for amplitude, which 
further depends on the masses of the planets, which are 
not known at this initial stage of modeling. In contrast, 
the dominant frequency from those models should be 
present in the data, although perhaps no longer domi- 
nant, even if both planets are eccentric or of a mass dif- 
ferent than assumed. This test for a TTV periodicity at 
the predicted frequency may not yield a significant detec- 
tion if (a) eccentricity significantly changes the character 
of the variations, or (b) the transit variations are dom- 
inated by a non-transiting planet. In these cases, large 
variations may be present, but have little power at the 
predicted frequency, so more data-driven methods (pa- 
pers II and III) are more useful. However, we note that 
for the systems presented in those papers, the dominant 
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Figure 1. Computed and measured durations of tlie planet can- 
didates in our sample. The line shows 1:1. The data come close 
to it, providing a useful check on our interpretation of planetary 
systems orbiting the targets stars. That the data fall slightly below 
the line could mean either the stellar parameters or the assumed 
edge-on and circular orbits are not perfectly correct. 



timescale of variation has a simple physical interpreta- 
tion and is captured by the nominal models described in 
this paper. 

The amplitude of the best-fitting sinusoid at the the- 
orized frequency is recorded, and it is compared to the 
same calculation for 10^ realizations of the observational 
O — C data with scrambled entries. These scrambled 
datasets should not contain the signal, but they will con- 
tain the same distribution of noise (which is possibly 
non-Gaussian, with outliers). Our criterion for confir- 
mation is that fewer than 100 of these realizations yield 
larger amplitudes than the real data, i.e., a False Alarm 
Probability (FAP) of < 10"^. Note that our method 
does not require that both planets show a TTV signal 
with such a low FAP. Even if the transits of one planet 
are at the limit of detectability, as long as we can es- 
tablish its orbital period and phase, we have what we 
need to estimate its dynamical effect on the other planet, 
whose TTVs must behave as expected. The FAP statis- 
tics are reported in Table 2] for the systems confirmed 
principally via this method ( i.e., those of Tables [l]|3l) . 
More over, we present sm all FAPs recon firming Kepler- 
9bc (IHolman et al.l [20Tol ). Kepler-18cd ([Cochran et al.l 
120111 ). and nearly all of the systems presented in papers 
II and III, as these systems show anticorrelated TTVs at 
the timescale predicted by the nominal models. 

3. RESULTS 

3.1. Confirmation of Kepler Planets 

Here we describe the 4 planetary systems confirmed 
principally by this method. Where interesting compar- 
isons can be made to other methods (paper II and paper 
III), we do so. Table [T] gives the catalog properties and 
Table [2] gives the physical properties of the host stars 
detailed in this paper. 

3.1.1. KOI-738 = Kepler-29 
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The lightcurve of Kepler-29 is plotted in Figure [21 and 
the phase curves of the two Neptune-size planets is plot- 
ted in Figure |3l Their periods of 10.3376(2) days and 
13.2907(4) days form a ratio 1.28566±0.00005, which lies 
right at a 9:7 {— 1.28571) orbital resonance. Being frac- 
tionally closer than 8 x 10^^ of a second-order resonance 
in this region would happen only 0.5% of the time by 
chance, and similar for other resonant orders. Being so 
close to a low-order resonance — consistent with dynam- 
ical libration — suggests dynamical interaction helped to 
establish their orbital architecture. Thus, from the pe- 
riod ratio alone, we have an indication that these objects 
are in the same system. 

Furthermore, since the ratio is not significantly offset 
from the resonance, the nominal integrations result in 
an extremely long-timescale variation in O — C (fig. [3]), 
even slightly longer than the sample 10-year integration 
used to find it. Therefore, on the timescale of our data, 
the predicted sinusoid would look only like a quadratic in 
O — C, as is observed (figure 01 left-hand side). The data 
for the shorter-period planet (Kepler-29b) indicate this 
quite robustly {FAP ~ 0.02%) and hint at oppositely- 
directed curvature for the longer-period planet (Kepler- 
29c has FAP ~ 1.5%). On the basis of the former we 
confirm the system really is composed of two planets. 

However, a parabolic signature is not specific enough 
to uniquely specify the other planet as the perturber, so 
we carried the investigation one more step. The nominal 
prediction for year-timescale O — C (figure |4l right-hand 
side) is not dominated by a parabola, but is instead a 
"chopping" signal on the frequency at which the planets 
come back to the same configuration, 

Pt4^ = 9 X 10.3376d = 7 x 13.2907d = 93d. (1) 

Thus, we search for that timescale. After subtracting 
the curvature models shown in figure [H (left-hand side), 
we computed the amplitude of the best-fitting sinusoid at 
the period of equation[Tl For Kepler-29b and Kepler-29c, 
the amplitude was 5 and 12 minutes, respectively, quite 
similar to that predicted by figure 01 Then we scrambled 
the data to find a FAP: it was 0.6% and 1.1% for the two 
planets, respectively. To emphasize this point, we plot in 
Figure the cross-correlation statistic defined in Paper 
III: a spike is seen at ^ 90 days, at positive values of S 
(i.e., anticorrelation at the predicted timescale). There- 
fore, the signal at the expected chopping timescale is a 
third, albeit modestly significant, indication that these 
two objects are really planets in the same system. 

3.1.2. KOI-806 = Kepler-30 

Kepler-30 is a remarkable system of three planets, all 
with extremely significant transit timing variations. Its 
lightcurve is shown in Figure [Bland each planet's phase 
curve is shown in Figure [71 A Neptune-size inner planet 
(Pfe = 29.2 days) lies just interior to a 2:1 resonance with 
a Jupiter-size planet {Pc — 60.3 days). Their mutual 
perturbation causes a transit timing variation with a fre- 
quency equal to the difference of the orbital frequencies 
from the resonance: 

P^^^ = l/|2/60.3d- l/29.2d| = 900d. (2) 

Since one of the participating planets has deep transits 
and thus is presumably massive, the signal is expected to 



be large. The theory predicts a sinusoid with a period a 
bit longer than the data. Indeed, the data show an enor- 
mous signal consistent with these expectations (Fig. [8]). 
Thus, we confirm Kepler-30b and Kepler-30c as objects 
orbiting the same star. 

Since their orbits are spaced a bit wider than the other 
planets confirmed in this paper, dynamical stability is 
not as strong a constraint on their masses. We find a 
mass upper hmit (see section 14. ip of about 17 Mjup. 
However, eccentricity can also account for the large tran- 
sit timing signals, so neither planet is necessarily nearly 
this massive (in particular, the small-radius Kepler-30b 
should not be multiple Jupiter masses). The devia- 
tions from a constant ephemeris for planet b are roughly 
parabolic and span ±10 hr, by far the largest timing vari- 
ations for a multiple-planet system yet detected (but not 
quite rivaling the circumbinary planet Kepler-16(AB)b; 
iDovle et al]l2011h . However, there is more signal in this 
TTV curve than just the large parabolic trend. The ris- 
ing branch of both the observed signal and the theoret- 
ical signal show a short-timescale "chopping" effect for 
Kepler-30b, where the even and odd transits times show 
structure, due to the 2:1 resonance. A similar effect was 
seen previously in Kcpler-9b (jHolman et al.l [2010[ ). A 
mass limit that takes into account the observed transit 
times will be discussed in section It results in a much 
tighter constraint, but it is less robust because it is not 
based on an exhaustive exploration of parameter space. 

We can also confirm Kepler-30d is likely a planet in the 
same system, as transit timing variations with a '-^ 2 hr 
timescale are seen in it (fig.[n]). However, the interpreta- 
tion is slightly less clear, as only 4 transits of Kepler- 
30d (P = 143.2128 d) have been recorded so far, so 
the shape of the variation is poorly determined. This 
perturbation is anticorrelated with the swing seen in its 
neighbor (Kcpler-30c) , but that swing might be predomi- 
nantly due to Kepler-30c's near-resonant interaction with 
Kepler-30b. According to the nominal integration, the 
timescale for the interaction between planets d and c 
should be ^ 400 days (Fig. [9l and Table [t]), similar to the 
span of the dataset. When the dataset is doubled, this 
might be discerned from the longer timescale (eq. [2]) in- 
teraction between the inner two planets. Even without a 
detailed analysis, we determine from transit timing that 
Kepler-30d is an object in the same system. We carry out 
preliminary transit-timing modeling of the three planets 
in section [4?2l 

Due to the large timing variations and deep transits 
found in this system, we find it instructive to show 
individual transits to make several other points. In 
Figure [TOl we plot the transits of Kepler-30b. Each 
portion of the lightcurve is centered where the tran- 
sit would be if it followed the best-fit linear ephemeris. 
Instead, we see that the timing of the planet varies 
by of order a day . The standard Kepler pipeline 
([Jenkins et al.|[20Tg ) assumes near-periodic transits, but 
the transits are deep enough for detection despite its 
strong acceleration. If the transits were considerably 
shallower, it may not have yielded a significant signal 
before the phase changed by more than a tr ansit du- 
ration ([Garcfa-Melendo fc L6pez-Moralesll2011[ ). To ad- 
dress this, the transit-timing subteam of Kepler is cur- 
rently developing search algorithms that assume either a 
quadratic, a sinusoidal, or a non-specified quasi-periodic 
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Figure 2. Kepler-29 lightcurves. Upper panel: the quarter-normalized, calibrated Kepler photometry (PA); lower panel: the detrended, 
normalized flux. The transit times of each planet are indicated by dots at the bottom of each panel. 
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Figure 3. Kepler-29 phase curves based on detrended, normalized 
flux. Each transit is shifted to its measured midtime, and transits 
with midtimes within 5 hours of a different planet's midtime are 
excluded from both this plot and from the model fits. Overplotted 
are transit models, box-car smoothed to the 29.4 minute cadence. 
The colors correspond to the dots of figure [2] 



ephemeris. 

Finally, consider Figure [Til the individual transits of 
Kepler-30c. The residuals (with a 3x zoom) are shown 
beneath each transit. The average transit lightcurve 
does not capture some of the variations from transit-to- 
transit, at the level of ^ 10% of the planet's depth. We 
suggest this is due to the planet transiting over starspots, 
which are clearly visible in the out-of-transit curve at the 
level of a few percent (fig.lHl)- To extract the spot period 
from the lightcurve, we used a slightl y modified version 
of the Discrete Correlation Function of E delson fc KrolikI 
(|T988) (see also White & Peterson 199|), as was re- 
cently applied to the CoRoT-7 data (jOueloz et al. 2009), 
and obtained Pjot = 15.25 ± 0.25 days. Modeling the 
spot-crossing signal during transits may allow the mea- 
surement of the spin orientation of the host star rel- 
ative to the planets' orbits in this particularly valu- 



able multiple-p lanet s ystem (iSanchis-Oieda et all 120111 : 
iNutzman et"al1 120111: iDesert et all 120 111) . Kepler-30c 
(also known as KOI-806.02) was recently detec ted in 
transit from the ground, bv iTinglev et al.l (|2011l) . who 

made the point that it will serve as an interesting sys- 
tem to continue following even after Kepler finishes its 
mission. 

3.1.3. KOI-935 = Kepler-31 

Kepler-31 is a system of four planet candidates spaced 
just wide of a 1:2:4:8 resonant configuration. Its 
lightcurve is shown in Figure [T^] and its phase curves 
in Figure [T51 From the nominal models, only the middle 
two (Kepler-31b and Kepler-31c) are expected to show 
a > 1 min variation, big enough to be detected in Ke- 
pler data. The timescale due to this near-resonance is 
predicted as 

P^^.^ = l/|2/42.6330d- l/20.8609d| = 980d. (3) 

Thus, we have the now-familiar situation that a long- 
timcscale variation is predicted due to the interaction. 
In this case we supplemented the dataset with quarters 
7 and 8 as well, which were able to capture more of this 
long-timescale variation. The data for Kepler-31b indeed 
clearly show a variation matching the prediction (Fig. 1141 
top panels), with a low false alarm probability. We are 
seeing its orbital period increase due to the torque ex- 
erted by a 2:1 mean-motion resonance. The magnitude 
of this torque depends on both the planets' masses and 
eccentricities, so it is not surprising the theory underpre- 
dicted the amplitude by a factor of ^ 3. Thus we confirm 
that these two objects are interacting, which combined 
with the mass limits set in section |4T1 confirms them as 
planets. 

The hint (FAP=3.4%) of TTVs of Kepler-31c due to 
Kepler-31b does not itself satisfy our statistical criterion. 
Nevertheless, the clear detection of TTVs in Kepler-31b 
consistent with the timescale predicted due to Kepler- 
31c provides evidence that they are in the same physical 
system. 
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Figure 4. Timing diagrams of Kepler-29. In the left hand panels, the data are plotted. In red is the best-fitting sinusoid of the 
theoretically-determined dominant perturbation frequency. From the amplitude of this sinusoid we determine the significance of the TTV 
detection. In the middle panels a nominal model, in which the planets are presumed to start out on circular, coplanar orbits, shown for 
the decade-span on which we determine the dominant perturbation frequency. Although this is a poor fit to the data, it illustrates that a 
very long-timescale and large-amplitude variation is expected due to the planets' orbits being right on top of a resonance with each other. 
Also, a very much smaller variation with a 93 day period is expected. In the right-hand panels, we have zoomed in on the first part of that 
signal; these variations are marginally detected in the data. 



15 - 



- 




100 200 300 400 500 600 



Period (days) 

Figure 5. Cross-correlation statistic (paper III) as a function 
of TTV period for Kepler-29 b/c. The dominant variation is of 
long period (as in Figure |4]l, and a secondary and only marginally 
significant peak occurs at ~ 90days. 



The system also has a long-period candidate KOI 
935.03 and an inner candidate KOI 935.04, but these 
are not visibly interacting. 

3.1.4. KOI-952 = Kepler-32 

Kepler-32 boasts five planet candidates, one of th e rich- 
est sys tems among the Kepler discoveries ( Lissaue r et al.l 
I2011bf) . Its lightcurve and phase curves are plotted in 
Figures [13 and [TH respectively. The lightcurve has vari- 
ations suggesting a starspot signal, in which case we de- 
termine Plot = 37.8 ± 1.2 days using an autocorrelation 
technique as above. 



Candidates .01 and .02 have short orbital periods and 
lie close to a 3:2 mean-motion resonance: P.02/P.Q1 — 
1.483. Therefore, we expect a short TTV timescale: 

Pjt„ = 1/13/8. 7522d- 2/5. 9012d| = 260d. (4) 

The observed timescale of TTVs is consistent with 
this expectation from the nominal N-body integrations 
(Fig. I17p . The amplitude is somewhat larger than pre- 
dicted by the nominal model (~ 2.6 and 3.3 x). This 
is not unexpected as even modest eccentricities can sig- 
nificantly increase the amplitude and completely change 
the phase of the predicted TTVs. The amplitudes also 
scale with the planet masses, which could be larger than 
those adopted in our nominal model. The ratio of TTV 
amplitudes for Kepler-32b and Kepler-32c agrees with 
the predictions of the nominal model to ^ 20%. De- 
spite these uncertainties, the requirement of dynamical 
stability provides firm upper limits on the planet masses. 
Thus, we confirm this pair of planets and rename them 
as Kepler-32b and Kepler-32c. 

Based on the nominal N-body integrations of the five 
planets, we expect shorter timescales and smaller ampli- 
tudes for transit variations of KOI-952. 03, .04, and .05. 
Indeed, no significant transit time variations were seen in 
them, nor do they induce detectable variation in Kepler- 
32b or Kepler-32c, so we cannot confirm them via TTVs 
at this time. 

3.2. Identification of the Host Star 

In this section we discuss the probability that the plan- 
ets are actually hosted by a star other than the target. 
The area on the sky in which that other star is allowed is 
limited via lack of movement of the centroid during tran- 
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Figure 6. Kepler-30 / KOI-806 lightcurves. Upper panel: quarter-normalized calibrated lightcurves; the stellar rotation period of 
15.25 ± 0.25 days is manifested as out-of-transit variations caused by spots. Bottom panel: Detrended, normalized flux. The transit times 
of each planet are indicated by dots at the bottom of each panel. 
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Figure 7. Kepler-30 / KOI-806 phase curves, in the style of figure 



sits (|Torres at al.ll20lil Bryson et al. in prep). Because 
of the low space density of transiting planetary systems, 
it is generally rather unlikely (quantified for each candi- 
date below) that it is hosted by an unrelated background 
or foreground star, which is a minority contribution to 
the light of the target. In comparison, because most stars 
are components of binaries, it is not inconceivable that 
our targets are physically bound binary stars where the 
second component is not detected, and the fainter com- 
ponent hosts the interacting planets (probability quan- 



tified below). To pursue these calculations, we assume 
a maximum planetary radius consisting of the upper en- 
velope of the mass-radius relation of known exoplanets. 
If the radii implied by a blend scenario are either above 
that limit or imply masses too large for dynamical stabil- 
ity given the configuration of planets, then we reject that 
blend scenario. Now we discuss details of each system. 

3.2.1. Kepler-29 

Kepler-29 has low contamination from surrounding 
stars. The nearest stars to Kepler-29, listed in the Ke- 
pler Input Catalog (KIC), are KID 10358756, & Kp = 
17.5 mag star located 6.5 arcsec to the south, and KID 
10358742 located 9 arcsec away. We examined the 
UKIRT J-band image and obtained a new image from 
Faulkes Telescope North in the Sloan-r band (figure [T8| : 
neither showed evidence of additional companions to a 
magnitude difference of ~ 5 in to ~ 1 arcsec. We there- 
fore adopt the aperture contamination calculated using 
stars in the KIC alone. 

The planets pass the centroid test: the photocenter 
out-of-transit and during transit are consistent to within 
(2(7, la) for planets (.01, .02) respectively. This rules out 
out background stars as the planet host if they are further 
than Rc = (0.7, 0.8) arcsec from the target, where Rc is 
the 3-cr radius of confusion. Thus none of the background 
stars in the image can be the host. Next we computed 
the space density of background stars, to determine the 
chance that the planetary system is actually around a 
star other than the target. That is, we assume that the 
TTV signature robustly indicates this system is plane- 
tary, but we wish to see the probability it is around a 
star other than the one we have characterized. 

Stability limits (section 14. ip require that the planets 
cannot be too massive relative to their host, otherwise 
they would be unstable. In this case, the mass ratio of 
both planets to their host must be < 2 x 10"''. With 
such low masses, we consider li?jup to be a robust upper 
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Figure 8. Kcpler-30. Transit timing diagram with theoretical model compared, for the interaction between Kepler-30b and Kepler-30c. 
Panels are the same as in Figure [4] 
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Figure 9. Kepler-30. Transit timing diagram with theoretical model compared, for the interaction between Kepler-30c (top) and Kepler- 
30d (bottom). (Ignore the labels above the plots for now, as they will be finalized with a Kepler number; these plots are for c and d.) 
Panels are the same as in Figure [4] 
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yet reported in multiple-planet systems. 



limit for the radii of the planets. This sets the constraint 
on the maximum magnitude difference that the blending 
star can have. Any star fainter than this limit would 
not contribute enough to the combined flux to generate 
an eclipse as deep as observed (830 ppm). Based on the 
U-shaped light-curve, we assumed that grazing transits 
were excluded, such that (Rp/R-i,)^ is a good approxi- 
mation of the (undiluted) depth. We used Yonsei-Yale 
isochrones (3 Gyr, solar metallicity) to relate the maxi- 
mum allowed magnitude difference (A Kp) to the spec- 
tral type and mass of the blending object. 

The background host hypothesis splits into two distinct 
cases: the planets are larger and (1) orbit a physically 
unassociated star, or (2) orbit an unseen binary compan- 
ion of the target. 

For t he former, we used the Besancon model of the 
galaxy ([Robin et al.l 120031 ) to estimate the number of 
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Figure 11. Indi vidual transits of K epler-30c (KOI-806.02, the 
planet observed bv lTinglev et ani2011l'l. in the style of figure [Tol 
Below each transit is shown the residuals from the model, scaled 
up by a factor of 3 to better show the deviations. The largest vari- 
ations appear in eclipse and are correlated point-to-point. Large 
starspots are visible in the lightcurve; we hypothesize this excess 
variability comes from the planet transiting over them. 



background stars within the region allowed by the cen- 
troid shift result and the region in the allowed magnitude 
difference (or stellar mass) . The probability of having a 
background star in this region is 0.6%. If we assume 
(conservatively) that the more massive planets on the 
background star are as likely, a priori, as less massive 
planets on the target star, then the odds ratio that the 
target star hosts the planets rather than a background 
star is 170 : 1. 

Next we calculate the probability of the latter case, 
that the planetary system is orbiting an unseen stellar 
companion to the target star. The blends involving stars 
less massive than a 0.6Mq star would need large planets, 
which in the range of known exoplanets would be too 
massive to be stable. In the range O.6-O.8M0, we can 
use the (Sloan) r - (2 MASS) K color of the target to 
rule out such blends, as they would not correspond to 
the observed color. A blended star more massive than 
O.SMq is allowed, and we call this scenario a 'twin star' 
blend. It is possible to put limits on such blends because 
we do not see two stars in the images, nor do we see 
two sets of lines in the spectrum. The only remaining 
case is that the two twins have a large separation along 
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Figure 13. Kepler-31 phase curves, in tiie style of figure |3] For 
the small inner candidate KOI-935.04, the phase is with respect to 
a linear ephemeris, the data in that panel are binned together in 
phase. The vertical scale of that panel is 20% of the other panels. 



the line of sight (a small velocity difference) but not in 
the plane of the sky. If there were such a case, it would 
correspond to twice the light, so the undiluted depth of 
transits would be twice as large, so the planets would be 
^ larger. 

Even without emp lo ying that constraint , we use 
iRaghavan eFall (|2010[ ): iDuquennov fc Mavod (|1991f ) to 
determine the distribution (frequency /mass ratio) of bi- 
naries and find the probability of having an unseen com- 
panion star in the allowed mass range (0.8-1.0Mq) is 7%. 

For this target, which has a clean centroid and closely- 



packed planets that would go unstable if they were gas- 
giant masses, we find that background hosts are much 
less likely than the target star being the host. 

3.2.2. Kepler-SO 

We obtained an I-band image of Kepler-30 from Lick 1- 
m telescope, and found no stars besides those in the KIC, 
the closest of which is Kp = 19.8 mag KIC 3832477, 
7.6 arcsec to the ENE. A star ~ 3 arcsec to the E is 
marginally detected in J-band with UKIRT, suggesting 
it is > 5 magnitudes fainter than the target in Kp. It can- 
not host the 2% deep transit of Kepler-32c, and we ne- 
glect its additional contribution to the dilution. Further- 
more, since planets c and d produce deep transits, their 
centroid information is particularly valuable. Their tran- 
sits have a consistent centroid with the target to ~ Icr 
in either case, and they limit the distance of a putative 
blend hosting the system to be within Rc = 0.2 arcsec of 
the target. 

As a routine part of vetting planet candidates, the 
depth for odd numbered transits was compared to the 
depth for even numbered transits. For the other candi- 
dates reported in this paper, there is agreement to ~ 3cr; 
however in this system, they disagree by 16(t for c and 
4cr for d. In some cases, that would point to a near- twin 
blended eclipsing binary causing the events. In this case, 
we have already identified via figure [TT] that transiting 
over starspots is the probable cause of these depth vari- 
ations. 

Using the same procedure as for Kepler-29, we find the 
probability that an unassociated background star hosts 
the planets is ~ 2 x 10"'*. In this case, a physical binary 
companion cannot host the planets, as then their depths 
would be large despite dilution, and their inferred radii 
would be larger than any planetary radii thus far mea- 
sured. 

3.2.3. Kepler-31 

Kepler-31 has low contamination from surrounding 
stars; no stars are seen within 8 arcsec in a UKIRT J- 
band image, and the closest comparably bright star is 
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Figure 14. Kepler-31. Transit timing diagram with theoretical model compared. The data for planets b and c show curvature on a long 
timescale, as predicted by a theoretical model in which they torque each other due to the close 2:1 resonance. 
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Figure 15. Kepler-32 lightcurve. Upper panel: quarter-normalized calibrated lightcurves; the stellar rotation period of 37.8 ± 1.2 days 
is evident in quarters 5 and 6. Bottom panel: The bottom, middle and top rows of dots (blue, green, teal) are for Kepler-32b, Kepler-32c 
and KOI-952.03. (The short-period candidates KOI-952.04 and KOI-952.05 are not detected in individual transits, so no pointers are given 
for them.) 
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Figure 16. Kepler-31 phase curves, in the style of figure |3] For 
the small inner candidate KOI-952.05, the phase is with respect to 
a linear ephemeris, the data in that panel are binned together in 
phase. The vertical scale of that panel is 20% of the other panels. 



KIC 9347893, 9.4 arcsec to the west. Moreover, the cen- 
troid information has all transits coincident within la 
of the target. The transits cannot be hosted by a back- 
ground star further than Rc = (0.3,0.5,0.8) arcsec in 
the case of Kepler-31b, Kepler-31c, KOI-935.03 respec- 
tively. For KOI-935.04, the transits are too shallow for a 
constraining centroid analysis. 

Again pursuing probability calculations as above, the 
chance of a star unassociated with the target being the 
actual host is only ^ 3 x lO^''. The probability of a 
physical companion hosting the planets is ^ 0.04. 

3.2.4. Kepler-32 

A J-band image from UKIRT shows the nearest star to 
be KID 9787232, - 6.6" to the west, resuhing in rather 
low contamination. 

The centroids during transit for Kepler-32b and 
Kepler-32c differ from those out-of-transit by only ~ 2a, 
roughly consistent with measurement uncertainties. The 
~ 3a radii of confusion Rc are 0.5" for Kepler-32b and 
0.8" for Kepler-32c. For KOT952.03, .04, and .05, the 
transits are too shallow for a constraining centroid anal- 
ysis. 

The host star is an M-dwarf and therefore of special in- 



terest. The Kepler Follow-up Program has obtained two 
spectra of Kepler-32: one spectrum from McDonald Ob- 
servatory and one from Keck Observatory. Both spectra 
are weak due to the faintness of the star (Kp=15.8). The 
cross correlation function between the observed spectra 
and available models is maximized for temperatures of 
^ 3900 K and ^ 3600 K, respectively. However, the 
atmospheric parameters are not well determined, as the 
star is cooler than the library of atmosphere models avail- 
able. Both spectra are consistent with the KIC clas- 
sification as a cool dwarf {Tcs = 3911, \ogg = 4.64, 
[M/II]=0.172). We conservatively adopt these values of 
Tcff and log g with uncertainties of 2 00K and 0.3 dex an d 
a [M/H] of 0± 0.4 based on the KIC (|Brown et al.llMl . 
By comparing to the Yonsei-Yale isochrones, we derive 
values for the stellar mass (0.58 ± 0.05Mo) and radius 
(0.53 ± O.O4i?0) that are slightly larger than those from 
the KIC. We estimate a luminosity of 0.06 ± 0.02 Lq 
a nd an age of < 9Gyr. 

iMuirhead et al.l ()2011[ ) have also obtained high- 
resolution IR spectrum of Kepler-32=KOI-952, finding 
a stehar T^s = 3726^^?, [Fe/H]= 0.04t°;°o- Int erpret- 
ing their data via Padova models (jGirardi et al.|[2002il . 
they inferred a considerably less massive and smaller star. 
We encourage further detailed analyses of the host star 
properties, as these have considerable uncertainties that 
directly affect the sizes and masses for the planets. 

The probability of a star unassociated with the target 
being the actual host is only ^ 3 x lO^'^. The probabil- 
ity of a physical companion hosting the planets is ~ 0.34. 
This latter number is relatively large in this case because 
all the transit depths are small, so they could in principle 
be much larger planets hosted by a star which is dramat- 
ically diluted. This opens up the possibilities for a very 
large range of companions (down to masses as low as 
~ O.IMq) that could host one or more of these objects, 
as long as transits near apocenter are invoked to match 
the durations (fig. [J). 

4. PLANETARY MASS LIMITS 

4.1. Dynamical Stability Analysis 

Many of the systems in this paper and its compan- 
ions (Papers II and III) are not completely solvable 
with present data; e.g., the gravitational interactions 
of the component planets do not yield unique solutions 
for their masses. Rather, there exists degeneracy be- 
tween the masses and eccentricities, as was the case for 
Kepler-9 before radial velocity constraints were applied 
(|Holman et al.n 2010). However, we constrain them to 
be in the planetary regime because the pairs of plan- 
ets all have small period ratios. In two-planet systems, 
a sharp boundary exist s between provably stable orbits 
(jMarchal fc Bozislll982f ) and orbits that are allowed to 
cross, according to energy and angular momentum con- 
servation. This boundary is when the separation of the 
planetary semi- major axes, aout — o-im exceeds a certain 
number (2-\/3 « 3.46, for coplanar, circular orbits) of 
mutual Hill spheres. 



(^in H~ ^out 

2 V 



(5) 



When the separation is only slightly closer than this, 
numerical integrations generally show the planets chaoti- 
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Figure 17. Comparison of measured transit times (left) and transit times predicted by the nominal model (right) for a system containing 
only Kepler-32b (top) and Kepler-32c (bottom) . Details are described in the caption to Fig. [4] 




Figure 18. FTN optical (Sloan r) image of Kepler-29. The scale 
is 1 arcmin on each side, and north is up. 

cally interact and Irave close encounters after only several 
million orbits, an astrophysically short ani ount of time 
(|Gladmanl 119931 : iBarnes fc Greenberd I2006D . An excep- 
tion is that planets locked in mean motion resonances 
may avoid each other due to the correlated phases of 
their orbits; Neptune and Pluto are a familiar example 
(jCohen fc HubbardI [19651 ). 

We ran a suite of numerical integrations for each 
planetary pair with a significant TTV interaction to 
determine when the masses would be too large to 
yield a long-term stable system. We use d the HY- 
BRID integrator within the Mercury package (jChamberd 
[1999). Since we stopped after close encounters, the 
mixe d-variable symplectic alg orithm was exclusively 
used (iWisdom fc Holmanlll991l) . with a timestep l/20th 
of the inner planet's orbital period. An implemen- 
tation of general relativistic precession was included 



(jNobih fc Roxburgh! [19861: ILissauer et all I2011b[ ). We 
picked the planets' ratio of masses from their ra- 
tio of radii ac cording to the relat ion Min/Mout — 
{R,nlRoutf -°^ fLissaue r et aI1l2011bf) . We performed 6 
integrations with differing total mass, starting at the the- 
oretical stability boundary and becoming more massive 
by logarithmic steps of 0.25 dex. The orbits were taken 
as initially circular and coplanar, with phases determined 
from the Kepler data. 

We stopped the integration when planets came within 
3 ru of each other, calling this time the instability 
timescale. We ran only one integration at each mass, 
so we recognize that this investigation only identifies the 
order of magnitude of that timescale - however, consis- 
tent with previous invest igators (Ch ambers et al.l[l996l : 
iDuncan fc Lissaueil 119971 ) we find that a factor-of-two 
change in mass causes many orders of magnitude of dif- 
ference in instability time. Figure [T^ gives the instability 
timescales as a function of planetary mass. We determine 
the mass upper limits based on the minimum mass for 
which the integrations start going unstable on < 10 Myr 
timescales; for this purpose, we do not need to run the in- 
tegrations for the likely lifetimes of these systems. These 
limits are given in Table [3] These are quite conserva- 
tive upper limits, and the true masses could all be quite 
smaller. Full modeling of the systems, including eccen- 
tricities, is required for a true mass measurement via 
transit timing variations (e.g. Kepler-11 by Lissauer et 
al.). Nevertheless, this exercise shows that all of the sys- 
tems we are presenting should be considered "planetary" 
systems, rather than stellar systems. 

In this investigation, we have simulated only pairs of 
planets and neglected additional planets, either other 
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Figure 19. The time until instability, as a function of planetary 
mass. This figure shows that for these objects to be stable, they 
must have masses traditionally associated with the planetary do- 



transiting candidates or completely unknown planets. 
When third planets are added, stability constraints have 
been found by numerical integrations to become even 
more stringen t (iChambers e t aL 1996 ; iChatteriee et al.l 
I2008t [Fa brvckv fc Murrav-Clav 2010al). An interesting 
exception~was discussed by iRavmond fc Barnes! I2005L 
but the third planet in that case quenched secular ec- 
centricity cycles, which is not the chaotic mechanism for 
instability we have investigated above. Therefore the fu- 
ture confirmation or discovery of such planets will not 
compromise the conclusions here. 

Similarly, we neglected possible eccentricity in the 
planets. For non-resonant cases, non-zero eccentricities 
would only serve to bring the planets closer tog ether at 
conjunctions, making th em more unstable (Dunca n et al.l 
119891 : iZhou et al1 l2007^. However, eccentricities can ac- 
tually increase the stability of resonant pairs. This effect 
could prolong the lifetime of Kepler-29, such that its true 
upper-limit masses are underestimated by at least a fac- 
tor of several. Since these orbits are closely packed how- 
ever, resonances besides the 9:7 can overlap with it, lead- 
ing to chaotic insta bilities. Using the Aa/a > 1.3/i^/^ cri- 
terion for stability ()Wisdomlll98ClD , the maximum planet- 
to-star mass ratio is ^ = 0.02. If that mass is shared 
among the planets, then they both fall in the planetary 
regime (< lOMjup). 



4.2. Preliminary Dynamical Fits to the Transits Times 

We have used the method first develo ped for the 
Kepler-9 and Kepler- 11 dis covery papers (H olman et al.l 
I2OIOI : iLissauer et al.l[2bllal ) to fit preliminary dynamical 
models to all the planetary pairs of this paper and the 
two companion papers. 

Using the Levenberg-Marquardt algorithm to drive 3- 
body numerical integrations, minimizing the of the 
residuals of the data minus the model. The free param- 
eters are the mass Mp of each planet and the Jacobian 
coordinates at dynamical epoch BJD 2455220.0 of each 
planet: orbital period P, epoch of mid-transit Tq, and 
the in-sky-plane and perpendicular-to-sky-plane compo- 
nents of the eccentricity vector, ecosw and esinw, re- 
spectively. We ignore the dynamical effect of any plan- 



ets that remain candidates in this work, just focusing on 
the interaction between the planets that show significant 
TTVs. 

The resulting fit parameters may be found in Table [5l 
In many cases, the eccentricities and masses are very 
highly correlated, resulting in poor errors on each quan- 
tity. However, we report many more decimal places be- 
yond what seems significant, as to make these fits re- 
producible by others. We note the contribution to for 
each planet and the number of degrees of freedom (d.o.f.) 
for that planet, meaning its number of data points minus 
the 5 free parameters used to model it. (In two-planet 
fits, there were 10 free parameters total.) In fits that are 
statistically acceptable, should equal d.o.f. to within 
a few times V d.o.f.. Of the systems confirmed in this 
paper, one of the poorer fits was for Kepler-30 / KOI- 
806, which might be attributable to starspots' effect on 
estimates of transit times. 

We found that Kepler-29, the system at the 9:7 res- 
onance, needed special treatment. Allowed to fit freely, 
the eccentricities solved for large values, causing the plan- 
ets to cross. Because of phase protection by the reso- 
nance, this was allowed by the data, but after a secular 
timescale the planets began chaotic scattering. There- 
fore in the fit reported in Table 5, we restricted the ab- 
solute value of each eccentricity components to be less 
than 0.05. The result was integrated using the Bulirsch- 
Stoer algorithm in Mercury (Chambers 199 9) for 30 Myr, 
during which it showed stable and regular orbital oscilla- 
tions. In other systems the eccentricities are quite mod- 
erate compared with the separation in semi-major axis, 
and we have not verified stability for these. In future 
detailed fits to the data, stability could be used as a 
principle guiding the results. 

In only a few cases are the masses meaningfully mea- 
sured, according to the formal error bars (e.g. > 3a). 
These error bars sample only the local minimum of the 
fit. We recognize that drawing meaning from the lo- 
cal curvature about the entire probability surface is haz- 
ardous. However, we hope this preliminary work on tran- 
sit fits will inspire other investigators to exhaustively 
explore orbital configurations that fit the data. In the 
meantime, we tried to set 3ct upper limits on the masses 
of the planets via the method developed in lHolman et al.l 
([2010): we moved the mass of one of the planets away 
from the best-fit, incrementing by 50%, and solving 
for the other parameters each time. When had grown 
to 9 greater than the of the best-fit, we used the mass 
of the planet in that fit to define the 3(7 upper limit. 
These limits are reported in table [5] In comparison to 
the stability study, these TTV limits on mass are much 
tighter. However, this method is rather delicate, in that 
we have not fully explored the global parameter space 
for a possibly more massive planet. The stability study 
provides more robust mass upper limits on the masses of 
these planets. 

5. DISCUSSION 
Taking stock of the results, we have 

• developed a new approach to confirming planets, 
by testing whether the transit timing signal is con- 
sistent with that produced by a known perturber, 
a second transiting planet. 
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• applied this approach to Kepler transit timing 
data, and confirmed 9 planets in 4 planetary sys- 
tems (and reconfirmed 16 planets in 8 additional 
systems), 

• showed that their masses must be in the planetary 
regime via stability arguments, and provided dy- 
namical fits to the data. 

The systems discussed herein have remarkable proper- 
ties. 

Kepler-29 is clearly engaged in a second-order res- 
onance, which has only been ob served for the 3:1 
resonance previously ( HP 60532; iDesort et al.l 120081 : 
iLaskar fc Correial I2009D . In such cases, adiabatic cap- 
ture from low-eccentricity is possible, but at a fi- 
nite speed of migration, the capture probability is en- 
hanced if pl anets migrated towards ea c h other on eccen- 
tric orbits (|Rein fc Papaloizoul 120101: iMustill fc Wvatti 
120111). Prev iously for systems of super-E arths, theorists 
(iTerauem fc Papaloizou 2007; McNeil fc NelsonI [20Tol: 
llda fc Linll2010l: ILiu et al.ll2011t iPierens et al.ll2011[ ) have 
suggested planets could be about this close to one an- 
other and trapped in resonances, but they have always 
focused on first-order resonances. Within the context of 
those formation models, the relative proportion of higher 
order resonances is now a pertinent question. 

Kepler-30 is perhaps the most dramatic system de- 
scribed here, having properties similar to multiple planet 
systems long-known from Doppler surveys. The inner 
planet b has a > 1 day swing in its transit times, per- 
turbed by its near-2:l outer compa nion, c. Thi s situa tion 
is very similar to that predicted bv lAgol et al.l ()2005D for 
the 2:l-resonant pair GJ 876 b/c, as the outer two planets 
clearly have gas-giant size. Continued monitoring of the 
transits may be able to give unique parameters for this 
system, and future information about duration variations 
or lack thereof may allow a measurement of their mutual 
orbital inclination. Planet c has a clear spot-crossing sig- 
nal; its depth varies as a function of the spot phase of 
its host star. With careful spot modeling, this may al- 
low the interpretation of the orientation of the host star's 
spin relative to the planets' orbits. Moreover, its tran- 
sits are rather deep, from which we infer a large radius 
(Table [3l 1.27 ± 0.16i?jup). This is larger than expected 
for a planet receiv ing so litt l e irra diation: the theoreti- 
cal models by Fort nev et al.l (l2007f) are typica lly smaller 
than ^ 1.2i?jup, and lDemorv fc Seagerl (|2011l ) confirmed 
that observationally using Kepler giant-planet candi- 
dates. We entertained the notion that this increase d 
depth might be due to rings (jSchlichting fc Chang||20ll] ). 
but unfortunately the characteristic "shoulders" were not 
seen in the transit lightcurve. 

Kepler-31 has a 20.9-day planet and a 42.6-day planet 
now confirmed. Two other candidates, at 9 and 88 
days, are present yet remain unconfirmed. If they 
can be confirmed at a later date, this system will be 
in an extraordinary near-l:2:4:8 resonant chain. Evi- 
dence for that sort of architecture has been building 

on several d ifferent front s: radial veloci ty detections 

(HP 40307: IMavor et al.l [200l GJ 876: IRivera et al.l 
[2OT0t HP 20794: iPepe et al.l l20Tll ). and a dramatic 
4-planet system discovere d by direct i maging (HR 
8799: iFaijrvckv fc Murrav-Clavl l2010bl: iMarois et ahl 



l2010f ). It has even been suggested that the Solar Sys- 
tem gia nt planets began in a multi-resonant config- 
uration (iMorbidelh et al.l l2007t iThommes eFall I2008D . 
albeit with links more compact than 2:1 reson ances 
(|Masset fc Snel lgrove'200ll: IPierens fc Nelsonl[200l . 

Kepler-32 is a particularly clean case of anticorrelated 
transit times that also have a timescale matching base- 
line expectations. It is particularly interesting to con- 
firm these planets, because their host star is an M-dwarf, 
around which small planets appear to be particularly 
abundant (Howard et al. 2011). The lightcurve boasts 
3 additional candidate planets as well, making a partic- 
ularly rich system. 

We can attribute most of these signals to a slight dis- 
placement of the planets' mean motions from a strict 
commensurability: they are involved in a "great inequal- 
ity," the orbital element oscillations of Jupiter and Sat- 
urn due to their displacement from the 5:2 resonance. 
The same sort of perturbations are d etected in the sys- 
tem of planets orbiting PSR1257+12 (iRasio etall 119921: 
IMalhotra et al.lll992l: [Peaia[l993l: lmTszczanlll994D . The 
planets torque each other's orbits, resulting in a period 
increase or decrease, depending on where the location of 
conjunctions is relative to their periapses and to the line- 
of-sight. The departure from commensurability causes 
this location to move, generating a predictable fluctu- 
ation in their orbital periods. The timescale for this 
fluctuation is easily calculated from their orbital peri- 
ods: equations [l]|31 A general feature is that large- 
amplitude timing signals take many orbits to manifest 
themselves. If the Kepler mission is extended to eight 
years, this method would reach its full potential for plan- 
etary pairs with longer "great inequality" or libration 
timescales (e.g., Kepler-29, flg.S]). As the Kepier mission 
seeks to confirm longer-period candidates, in particular 
candidates in the habitable zone, we will be attempting 
transit timing analyses based on fewer transits. Pynami- 
cal theory may be required to condition our expectations 
about transit timing variations. This paper is a stepping 
stone to that type of analysis. 

Along with Papers II and Paper HI, we have confirmed 
21 planets in 10 systems. The transit timing method has 
very little bias with respect to the magnitude of the stel- 
lar host, as planetary systems with large, detectable per- 
turbations are hosted by stars throughout the sample. 
This contrasts with other Kepler confirmations so far, 
which have mostly relied on ground-based (and Spitzer 
Space Telescope) follow-up to confirm the planetary na- 
ture of the transiting objects. In particular, we show in 
Figure that the current TTV confirmations are drawn 
from a very much fainter population than the previous 
confirmations. 

In two of these systems, as well as systems in Papers 
II and HI, additional candidate planets have been found, 
but not confirmed via TTV or other methods. A priori, 
it has been argued that planet candidates in multiple- 
planet systems have a higher fidelity than planet can- 
didates that are by themselves, e ven before performing 
follow-up observatio ns or analysis (jLissauer et al. 20111a; 
iLatham et al.l[20Tll Lissauer et al. in press). Previously 
our team has pursued further analysis to validate ad- 
ditional s mall planets in TTV -active plane tary systems: 
Kepl er-9d (jTorres et al.ll201lD. Kepler-lOc (iFressin et al.l 
[20ll . and Kepler-18b ICochran et aIl l201lD. For the 
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Figure 20. Histogram of stellar magnitudes in ttie Kepler band 
(Kp). The TTV catalog of planets confirmed by their interaction 
with each other (the 10 planetary systems from Paper II, Paper III, 
and this paper) is not dependent on follow-up with other telescopes, 
so it more uniformly samples the intrinsic magnitude distribution 
of Kepler stars. The previously confirmed pla nets, by t he Kep Jer 
team (Kepler - 4-Kep l er-22) and also by Santcrnc et al.l (I2011b|[3); 
IBouchv et aP |(20T]] '| : [Bonomo et al.l pOllI: Johnson et al. (2011|), 
do depend on other telescopes, and thus tend towards brighter 
targets. 



candidates listed here, many have good hmits on Rc, the 
distance to which an unseen blend is a possible source of 
these transits. Their phased lightcurves are easily fit by 
planetary lightcurves ("U" shaped) rather than prefer- 
ring blended binary lightcurves (which are usually "V" 
shaped). Finally, the additional candidates have dura- 
tions that can be explained by orbiting the same stars as 
the confirmed planets (Fig.[T]). Despite the fainter target 
stars which makes formal validation difficult, we expect 
most if not all of these candidates are indeed planets. 
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observations obtained with facilities of the Las Cumbres 
Observatory Global Telescope. 

Facilities: Kepler. 



Confirmation of 4 Multiple Planet Systems 



Table 1 

Properties of target stars 





KOI 


KIC-ID 


Kp 


CO" 


CI" 


C2" 


C3" 




RA 


DEC 
















[ppm] 


hr [J2000] 


deg [J2000] 


Kepler-29 


738 


10358759 


15.282 


0.108 


0.054 


0.091 


0.092 


176 


19 53 23.60 


+47 29 28.4 


Kepler-30 


806 


3832474 


15.403 


0.094 


0.099 


0.116 


0.056 


652 


19 01 08.07 


+38 56 50.2 


Kepler-31 


935 


9347899 


15.237 


0.053 


0.045 


0.119 


0.063 


186 


19 36 05.52 


+45 51 11.1 


Kepler-32 


952 


9787239 


15.913 


0.096 


0.117 


0.193 


0.119 


288 


19 51 22.18 


+46 34 27.4 



Note. — Information mostly from the Kepler Input Catalog. 
^ Contamination for each season 0-3 ( season — (quarter-|-2) mod 4 ): 
the fractional amount of light leaking in to the target's aperture from 
other stars, known from the Kepler Input Catalog. 



Table 2 

Table of Stellar Properties of Hosts 







logs 


vsini 


[Fe/H] 




R* 




(K) 


(cgs) 


(km s-l) 




Me 


Rq 


Kcplcr-29 


5750 ± 250 (L) 


5.00 + 0.25 (L) 


4 + 2 (L) 


0.0 + 0.3 (N) 


1.00+0.12 


0.96+0.14 


Kepler-30 


5498 ± 54 (K) 


4.77 + 0.23 (K) 


1.94 + 0.22 (K) 


0.18 + 0.27 (K) 


0.99+0.08 


0.95+0.12 


Kepler-31 


6340 ± 200 


4.696 ± 0.300 


-0.076 + 0.400 


1.21+0.17 


1.22+0.24 


Kepler-32 


3900 ± 200 (K,M) 


4.64 ± 0.30 




0+ 0.4 


0.58 + 0.05 


0.53 + 0.04 



^ Sources for stellar properties. Spectroscopic parameter with un- 
certainties indicated in parentheses are from: K— Keck Observatory, 
L— Lick Observatory, M^McDonald Observatory, N— NOAO. Quoted 
uncertainties do not include systematic uncertainties due to stellar 
models. 
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Table 3 

Key Properties of Planets and Planet Candidates 





To"" 


pb 




fed 


U d 


Depth 


Rp/ i^*^ 




a 






(d) 


(d) 


(d) 






(ppm) 


R^ 


(AU) 




29b 


82.750±0.009 


10.3376±0.0002 


0.117±0.003 


0.0 - 0.85 


0.0 - 0.5 


1204 


0.0343±0.0008 


3.6±0.5 


0.09 


0.4 


29e 


78.471±0.016 


13.2907±0.0004 


0.127±0.006 






871 


0.0280±0.0012 


2.9±0.4 


0.11 


0.3 


30b 


83.04±0.41 


29.329±0.022 


0.200±0.004 






1540 


0.0351±0.0005 


3.6±0.5 


0.18 


0.2 


30c 


176.904±0.005 


60.3251±0.0008 


0.2392±0.0011 


0.44±0.04 


0.584±0.023 


20578 


0.1375±0.0011 


14.3±1.8 


0.30 


9.1 


30d 


87.220±0.038 


143.213±0.013 


0.322±0.003 


0.52±0.04 


0.56±0.05 


11014 


0.1020±0.0014 


10.6±1.4 


0.5 


17 


935.04 


85.10±0.04 


9.6172±0.0005 


0.176±0.014 




- 0.6 


161 


0.0136±0.0006 


1.8±0.4 


0.09 




31b 


92.141±0.006 


20.8613±0.0002 


0.208±0.002 


- 0.75 


0.38 - 0.70 


1895 


0.0411±0.0006 


5.5±1.1 


0.16 




31c 


74.191±0.007 


42.6318±0.0005 


0.251±0.004 


- 0.8 


0.18 - 0.67 


1729 


0.0400±0.0007 


5.3±1.1 


0.26 


4.7 


935.03 


67.942±0.009 


87.6451±0.0014 


0.344±0.010 


- 0.9 


0.1 - 1.0 


959 


0.0291±0.0011 


3.9±0.8 


0.4 


6.8 


952.05 


65.54±0.04 


0.74296±0.00007 


0.039±0.004 






224 


0.0142±0.0007 


0.82±0.07 


0.013 




952.04 


66.61±0.03 


2.8960±0.0003 


0.053±0.004 






671 


0.0259±0.0013 


1.5±0.1 


0.033 




32b 


74.902±0.008 


5.90124±0.00010 


0.088±0.005 






1650 


0.0389±0.0019 


2.2±0.2 


0.05 


4.1 


32c 


77.378±0.013 


8.7522±0.0003 


0.097±0.010 






1453 


0.0352±0.0033 


2.0±0.2 


0.09 


0.5 


952.03 


88.211±0.011 


22.7802±0.0005 


0.124±0.003 




0.0 - 0.75 


2181 


0.0467±0.0011 


2.7±0.2 


0.13 





^ BJD-2454900. Rather than statistical error bars, we report the RMS 
of the measured transit timing deviations. 

^ Derived from the measured transit times. Error bar is the RMS of the 
measured transit timing deviations, divided by the number of transits 
that the data span. 

Duration of time that the center of the planet overlies the stellar 
disk. 

Values with it error bars are from formal fits. However, often we 
found b and u to be quite poorly constrained and/or degenerate, so we 
computed grids in u G [0, 1], and give either a 2— a (Ax^ — 4) range, 
or if this range covers the whole grid, no value at all. 

Takes into account contamination values from Table fl] 
^ Using stellar radii from Table [2] 

^ Based on assumption of dynamical stability and stellar mass from 
Table [2] 
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Table 4 

Table of Statistics for Pairs of Planets Candidates 



KOIin frequcncyin FAPin KOIout frequcncyout FAPout 



Kepler-29b 0.00026034 



0.00018 

< 10-^ 
0.70891 
0.00855 

< 10-^ 
0.53060 
0.02580 



Keplcr-29c 0.00026086 0.01513 



Kepler-30b 
Kepler-30b 
Kepler-30c 



0.00098529 
0.00698174 
0.00259824 



Kepler-30c 
Kepler-30d 
Kepler-30d 



0.00099742 
0.00076809 
0.00262702 



0.00016 

0.29269 
< 10-^ 



Kepler-31b 
Kepler-31b 
Kepler-31c 



0.00100550 
0.01140540 
0.00064850 



952.04 
952.04 
952.04 
Kepler-32b 
Kepler-32b 
Kepler-32c 



Kepler-31c 
935.03 
935.03 

Kepier-32b 

Kepler-32c 
952.03 

Kepler-32c 
952.03 
952.03 



0.00100732 
0.00229586 
0.00064007 



0.03360 
0.96650 
0.13400 



0.00642701 
0.11678254 
0.04389705 
0.00390700 
0.04389928 
0.02646319 



0.11427 
0.09436 
0.51453 
0.00001 

0.84870 
0.87862 



0.00642994 
0.00252686 
0.00587133 
0.00391415 
0.00613203 
0.01743588 



Kepler-9b 0.00066337" 



0.00002 
■5— 



Kepler-9c 0.00066211 



0.27161 
0.16165 
0.61985 
0.05363 
0.49878 
0.65570 



0.00170 
< 10^ 



Kepler-18c 0.00373946 < 10 



168.02 
168.02 
Kepler-23b 



Kcplcr-18d 
Kepicr-23b 
Kepler-23c 
Kepler-23c 



0.00381809 



0.02919133 
0.01021819 
0.00215313 



0.80572 
0.37630 
0.00084 



0.02919325 
0.01022921 
0.00214191 



0.80219 
0.65446 
< 10"^ 



1102.04 


0.00994495 


0.22250 


1102.02 


0.00994282 


0.35748 


1102.04 


0.07345063 


0.75543 


1102.01 


0.00762697 


0.91916 


1102.04 


0.05264294 


0.25841 


1102.03 


0.02506453 


0.56294 


Kepler-24b 


0.00237619 


< 10"^ 


Kepler-24c 


0.00237709 


< 10"^ 


Kepler-24b 


0.01749593 


0.85933 


1102.03 


0.01752736 


0.14017 


Kepler- 24c 


0.00424632 


0.63164 


1102.03 


0.00424409 


0.57523 


Keplcr-25b 


0.00305506 


< 10" 


Kepler-25c 


0.00306491 


0.02136 


250.03 


0.11935114 


0.02223 


Kepler-26b 


0.03794286 


0.84931 



250.03 0.05797149 0.06283 

Kepler-26b 0.01119109 0.05656 

Keplcr-27b 0.00137000 O.OOOW 

Keplcr-28b 0.00426232 < 10-^ 



Kepler-26c 
Kepler-26c 
Kepicr-27c 



0.00765104 
0.01118508 



0.10796 
0.28181 



0.00137385 0.00083 



Kepler-28c 0.00429017 



0.00239 



Note. — Theoretically predieted O-C frequencies are in cycles 
per day. False Alarm Probabilities (FAP) are bold, if the detection 
is considered significant [FAP < 10^^). Systems above the line are 
presented as planetary system discoveries in this paper. Systems below 
t he double horizontal l ine are presented as plane tary system discoveries 
in lHolman et all J2010I ) , in lCochran et all 1 120111) , and in the companion 
papers, papers II and III. 
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Table 5 

Example TTV Solution for Planets Candidates 
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Note. — 

here. 



A sample of TTV-fit results for each system presented 
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Table 6 

Transit Times for Kepler Transiting Planet Candidates 



KOI 


n 




TTV„ 








BJD- 2454900 


(d) 


(d) 


738.01 




82.749505 -i 


- n X 10.337583 




738.01 





82.7642 


0.0147 


0.0081 


738.01 


1 


93.1033 


0.0162 


0.0074 



Note. — Tablc[6]is published in its entirety in the electronic edition 
of the Astrophysical Journal or from the "tttable" within the source 
files of this arxiv posting. A portion is shown here for guidance regard- 
ing its form and content. 
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